Equation of state for warm dense lithium: A first principles investigation
Long Feiyun1, Liu Haitao1, Li Dafang1, Yan Jun1, 2, †
Institute of Applied Physics and Computational Mathematics, Beijing 100088, China
Center for Applied Physics and Technology, Peking University, Beijing 100871, China

 

† Corresponding author. E-mail: yan_jun@iapcm.ac.cn

Abstract

The quantum molecular dynamics based on the density functional theory has been adopted to simulate the equation of state for the shock compressed lithium. In contrary to some earlier experimental measurement and theoretical simulation, there is not any evidence of the ‘kink’ in the Hugoniot curve in our accurate simulation. Throughout the shock compression process, only a simple solid-to-liquid melting behavior is demonstrated, instead of complicated solid–solid phase transitions. Moreover, the x-ray absorption near-edge spectroscopy has been predicted as a feasible way to diagnose the structural evolution of warm dense lithium in this density region.

1. Introduction

Lithium, as the first alkali metal element in the periodic table, seems to be very simple formally. For example, the lithium atom has only occupied s-shells, and the body-centered crystal structure is found as the ground state of bulk lithium at room temperature. However, nontrivial physical phenomena have been observed under the change of temperature and pressure, to attract long-standing research interests until now. For instance, abundant phase transitions are extensively studied and plenty of new phase have been found,[111] as well as the anomalous melting behavior has been also investigated.[1218] Moreover, to change the external conditions, the electric conductivity of lithium could be altered greatly, e.g., metal, non-metal, and even superconductor.[10,1929] Furthermore, not only the fundamental curiosity, but also plenty of applications do require an accurate investigation for the lithium in the regime of higher temperature and pressure, i.e., the so-called warm dense matter (WDM). Especially, advanced experimental techniques[30,31] and theoretical methods[3234] will give a boost to the knowledge of warm dense lithium.

Despite great efforts spared by former researchers, even the equation of state of shock compressed lithium is still not certain so far. Debate remains about the shape of the Hugoniot curve around the several kilo-kelvins and about 2–3 times compression relative to ambient density. In detail, the first experimental measurement on the shock compressed lithium was performed by Bakanova et al.[35] Generally, the pressure was observed to increase with the increasing compression along the Hugoniot curve. But this ascending tendency was significantly lowered in the region between about 1.2 and 1.4 for the density and a distinct ‘kink’ could be observed there. However, this experimental discovery was contradicted by a theoretical simulation performed by Young and Ross in 1984.[36] In fact, their theoretical results not only disapproved the existence of the so-called ‘kink’ in the Hugoniot curve, but also deviated from several sets of experimental results under larger compression. A noticeable matter is that the parameters in this empirical pseudo-potential calculation were fitted with the room temperature isotherms of lithium, so it is doubtful whether this theoretical method was accurate throughout the shock compression process. Several years ago, Kim et al. re-examined the shock Hugoniot curve with the electron force field (eFF) method[37] developed by themselves.[38] The ‘kink’ was then confirmed again, and they explained this phenomenon with two successive phase transitions (fcc-to-cI16 and cI16-to-amorphous). Their results are in good agreement with several groups of experimental data,[35,3941] though most of them were not available for the density larger than 1.2 . On the other hand, the phase transitions predicted by Kim et al. seems to be astonishing from the aspect of melting, because the melting point of compressed lithium is generally only several hundred kelvins.[14,15] In addition, the melting point is decreased with increasing pressure under large compression.

To clarify the long-standing confusion, it is necessary to investigate the equation of state of warm dense lithium with an accurate theoretical method independent with empirical parameters. In this article, we plan to adopt the quantum molecular dynamics (QMD) to overcome the inaccuracy from empirical parameters existed in previous theoretical simulations. Specifically, the shock Hugoniot curve is investigated at first, and then the ion–ion pair correlation function (PCF) is also calculated for each key point in the curve to reveal the structural evolution during the shock process. Furthermore, we tend to seek a common experimental measurement which could demonstrate the structural evolution clearly. Since the x-ray absorption near-edge spectroscopy (XANES) could be greatly influenced by the fast varied chemical environment in WDM, it is well accepted as a convenient experimental technique for WDM. Herein the Kubo–Greenwood formula based on the linear response theory is taken into account.

2. Computational method

The Vienna ab initio simulation package (VASP)[42,43] has been employed for the QMD simulation in this article. To simulate the movement of ions, the Born-Oppenheimer molecular dynamics is performed, and the Nöse–Hoover thermostat[44] is enabled to control the ion temperature in the canonical ensemble. A large supercell containing 128 Li atoms is chosen to mimic the real system as accurately as possible. In general, a molecular dynamics simulation should last for 8 ps–12 ps with the time step of 1 fs, and only the trajectory in the final 1 ps–2 ps will be exploited to calculate the EOS and other properties like the PCFs. On the other hand, to describe the movement of electrons, the finite temperature density functional theory (DFT)[45] is suitable for achieving a good compromise between computational efficiency and accuracy. The assumption of local thermodynamical equilibrium was imposed all along the shock Hugoniot, hence the electron temperature in Fermi–Dirac distribution was set equal with the ion temperature . To approximate the exchange and correlation functional, the Perdew–Burke–Ernzerhof (PBE) functional[46] based on the generalized gradient approximation is selected owing to its popularity in simulations of condensed matter and WDM. Moreover, the electron–ion interaction is described by the projector augmented wave (PAW) potential.[47] In our simulations, the PAW potential with three valence electrons is chosen for reasonable results under higher pressure. For such a large supercell, even the gamma point should be sufficient for the Brillouin zone sampling in the molecular dynamic simulation. In the electronic self-consistent calculation, one of the most important parameters is the cutoff energy which could decide the size and quality of the plane-wave basis set. Several different values have been examined and the 700 eV has been considered to obtain a good convergence for both total energy and pressure.

For the shocked lithium, the equation of state could be calculated with the Rankine–Hugoniot equations . Herein the subscripts 0 and 1 denote the initial (cold) and shocked states, respectively. The variables E, V, and P represent the internal energy, volume, and pressure respectively. The internal energy equals to the total energy in the DFT calculation plus the zero-point energy, and the pressure includes both the electronic and ionic portions. Our simulations started from the solid density (i.e., ) under the room temperature, and the body-centered crystal is used as the initial structure. Table 1 listed the principal Hugoniot points. Our simulation are in the range of 0.6 –1.5 for the density and from the room temperature to near 16 kilo-Kelvin for the temperature, so we believe that it could be meaningful in determining the existence of the above-mentioned ‘kink’ in the shock Hugoniot curve.

Table 1.

The density (ρ), pressure (P), and temperature (T) of the warm dense lithium along the principal Hogoniot. The particle velocity ( ) and shock velocity ( ) are also presented.

.

In order to predict the XANES spectrum, the impurity model[48] is taken into account by using the ABINIT code.[49] The comprehensive description about this computational method could be found in Refs. [50] and [51]. In short, one atom among the whole supercell is selected as the ‘absorber’ atom immersed into the surrounding lithium atoms in this scheme. For this ‘adsorber’ atom, a new PAW with a 1s core–hole is generated with the ATOMPAW code,[52] as well as the standard PAW still served for other atoms. It should be noted that the absolute edge position of the XANES spectrum could not be obtained owing to the limitation of this method. Moreover, ten snapshots are chosen from the dynamical trajectory in equilibrium to calculate the average XANES spectrum for each state in the shock Hugoniot curve.

3. Results and discussion

In Fig. 1, our simulated Hugoniot curve of the shock-compressed lithium is shown, and some earlier experimental[35,3941] and theoretical[36,38] results are also plotted therein for an extensive comparison. At the beginning of the shock compression process, most results are in good agreement, till the density was increased to 0.9 . There are several noticeable characteristics, i.e., the eFF simulation results are generally slightly higher than others, and the Bakanova’s experimental results tend to give smaller pressure in comparison with other experimental and theoretical results. In this region, our work is in excellent agreement with most experimental results.

Fig. 1. (color online) The shock Hugoniot curve of lithium.

In respect to the enhanced compression, the various equations of state tend to diverge from one another for densities above 0.9 . In general, Kim’s theoretical results agree more with some experimental results.[39,41] However, the theoretical work by Young and Ross[36] gradually predicts higher values than the others. Because they adopted the empirical pesudo-potential method with fitting parameters in accord with the room-temperature isotherm, it is plausible that their results behaved poorly when temperature and density became higher. In fact, the overestimation of the pressure could be understood as the consequence of such a stiff pseudo-potential. Our simulation results are lower than both these results for the density up to 1.2 . The fact should be mentioned that our results are still in agreement with Bakanova’s experimental results[35] in this region, and the experimental results given by the Lawrence Livermore National Laboratory[40] are also consistent with our results for the density up to about 1.0 .

Nextly, we should pay more attention to the region of higher density, i.e., around between 1.2 and 1.4 , because of the above-mentioned argumentation of the ‘kink’ in the Hugoniot curve. Unfortunately, there is still lack of abundant experimental results in this region to make an comprehensive comparison. The only available experiment in this region is performed by Bakanova et al.,[35] to confirm such a ‘kink’ in the Hugoniot curve as mentioned above. In the theoretical aspect, Young and Ross disproved the existence of the ‘kink’ and presented a smooth increase of the pressure with the enhanced compression in the region.[36] However, their results behaved poorly even for lower temperature and density, and it is doubtful whether their empirical pseudo-potential method could be applicable for this case. On the contrary, a newer simulation by Kim et al. still confirmed the existence of the ‘kink’.[38] Based on their developed eFF quantum electron dynamics, a much more distinctive ‘kink’ is found. In spite of the fact that this work agreed well with some experimental results, great uncertainty still exist because the eFF method bring in some fairly great approximations to facilitate the large-scale simulation of complicated systems, e.g., electrons are represented as Gaussian wave packets and the wave function is even treated as a Hartree product. For our results based on the accurate QMD, the Hugoniot curve is still very smooth in this region, to decline such a ‘kink’ again. Interestingly, our results agree with Bakanova’s experiment very well before the ‘kink’ appeared in these experimental results. Moreover, the curve given by our calculation is similar with that by the empirical pseudo-potential method,[36] though the latter tended to greatly overestimate the pressure for the same density.

To understand the origin of the difference in EOS, the internal structure of the warm dense lithium should be examined at first. The PCF, who corresponds to the average neighbors around a certain ion at a specific distance, could be a key issue to reflect the structural evolution during the shock compression process. For the convenience of comparing results of various density, the PCF is scaled by , and n is the number density.

As shown in Fig. 2, with the increasing temperature and pressure, the first maximum in the PCF tends to be broaden gradually. For states with lower temperature and density (i.e., and ), there are still remarkable peaks at larger distance, in corresponding to a long-range order. To increase both temperature and density, the peak at the scaled distance of about 4.5 is greatly broaden and then disappeared to imply a loss of long-range ordered structure. The peak at the scaled distance of about 3.5 is also decreased gradually. Thus a simple PCF only with a predominant peak around 1.7 is observed since the density became larger than 1.2 . Clearly, the evolution process is corresponding to a typical solid-to-liquid transition (i.e., melting) herein, in quite different with the eFF simulation results,[38] where a sudden change is shown and then explained with solid-to-solid phase transitions. Furthermore, another survey has been made for the cI16 initial structure with the density of 1.3 ). By performing a series of QMD simulation at several different temperatures, the PCFs of these equilibrium states are given in Fig. 3, respectively. For the room temperature, several peaks are observed, in corresponding to the crystal structure. Even for slightly higher temperature at about 500 K, a typical liquid behavior already appears, to show the fact that the melting is so easy for this density. Therefore, it seems that the solid-to-solid transitions predicted by Kim et al. are impossible at all.

Fig. 2. (color online) The PCFs of the warm dense lithium along the Hugoniot curve. The number on each plot is corresponding to the density in units .
Fig. 3. (color online) The PCFs of the warm dense lithium at a fixed density of 1.3 and under different temperatures. The cI16 crystal structure is adopted as the initial structure.

In recent years, the XANES has been well shown as a powerful tool to investigate the microstructure of compressed materials, e.g., heavy alkali metal.[53] As shown in Fig. 4, the simulated XANES spectra of the warm dense lithium along the Hugoniot curve are given. For the solid crystal structure under the room temperature, plenty of peaks could be observed. To increase the temperature and density gradually, the XANES spectrum tends to be more and more smooth. Some typical characters could be discovered for the K absorption edge. With the increasing temperature and density, the position of the K absorption edge should be shifted to the higher energy. In an early experimental measurement on lithium, the similar blue-shift according the increasing temperature has been also reported.[54] Moreover, the slope of the K absorption edge will become lower when the temperature and density are increased. The XANES spectrum is greatly sensitive to the variation of the warm dense lithium, then a high resolution XANES experiment could be very helpful to approve our simulation results.

Fig. 4. (color online) The simulated XANES spectra of lithium along the shock Hugoniot curve. It should be noted that the spectra are shifted to match the K absorption edge position of the initial structure with the experimental value.
4. Conclusions

In summary, we performed an accurate investigation on the shock compressed lithium, with the help of the QMD based on the finite temperature DFT, in an attempt to clarify the long-standing confusion about the equation of state. Our results presented a smooth Hugoniot curve for the density up to 1.5 , and then declined the existence of a ‘kink’ in the region between about 1.2 and 1.4 . Moreover, a simple melting behavior is also approved in corresponding to the smooth Hugoniot curve, instead of the complicated solid-to-solid phase transitions reported before. Actually, the warm dense matter is very complicated even for a simple element like Li, thus the first principles calculation could play an important role herein. Furthermore, the XANES is suggested as an efficient experimental technique for demonstrating the structural evolution of the warm dense lithium along the shock Hugoniot curve.

Reference
[1] Hanfland M Syassen K Christensen N E Novikov D L 2000 Nature 408 174
[2] Orlov A I Khvostantsev L G Gromnitskaya E L Stal'gorova O V 2001 J. Exp. Theor. Phys. 93 393
[3] Lindl J D Amendt P Berger R L Glendinning S G Glenzer S H Haan S W Kauffman R L Landen O L Suter L J 2004 Phys. Plasmas 11 339
[4] Maksimov E G Magnitskaya M V Fortov V E 2005 Phys.-Usp. 48 761
[5] Pickard C J Needs R J 2009 Phys. Rev. Lett. 102 146401
[6] Yao Y Tse J S Klug D D 2009 Phys. Rev. Lett. 102 115503
[7] Yang J J Tse J S Iitaka T 2010 J. Phys.: Condens. Matter 22 095503
[8] Kulkarni A Doll K Prasad D L V K Schoen J C Jansen M 2011 Phys. Rev. 84 172101
[9] Lv J Wang Y C Zhu L Ma Y M 2011 Phys. Rev. Lett. 106 015503
[10] Orlov A I Brazhkin V V 2013 JETP Lett. 97 270
[11] Naomov I I Hemley R J 2015 Phys. Rev. Lett. 114 156403
[12] Tamblyn I Raty J Y Bonev S A 2008 Phys. Rev. Lett. 101 075703
[13] Lepeshkin S V Magnitskaya M V Maksimov E G 2009 JETP Lett. 89 586
[14] Li D F Zhang P Yan J Liu H Y 2011 Europhys. Lett. 95 56004
[15] Rousseau B Xie Y Ma Y Bergara A 2011 Eur. Phys. J. B 8 1
[16] Guillaume C L Gregoryanz E Degtyareva O McMahon M I Hanfland M Evans S Guthrie M Sinogeikin S V Mao H K 2011 Nat. Phys. 7 211
[17] Schaeffer A M J Talmadge W B Temple S R Deemyad S 2012 Phys. Rev. Lett. 109 185702
[18] Arafin S Singh R N 2016 J. Phys. Chem. Solids 91 101
[19] Shimizu K Ishikawa H Takao D Yagi T Amaya K 2002 Nature 419 597
[20] Fortov V E Yakushev V V Kagan K L Lomonosov I V Maksimov E G Magnitskaya M V Postnov V I Yakusheva T I 2002 J. Phys.: Condens. Matter 14 10809
[21] Bastea M Bastea S 2002 Phys. Rev. 65 193104
[22] Kasinathan D Kunes J Lazicki A Rosner H Yoo C S Scalettar R T Pickett W E 2006 Phys. Rev. Lett. 96 047004
[23] Profeta G Franchini C Lathiotakis N N Floris A Sanna A Marques M A L Lüders M Massidda S Gross E K U Continenza A 2006 Phys. Rev. Lett. 96 047003
[24] Kietzmann A Redmer R Desjarlais M P Mattsson T R 2008 Phys. Rev. Lett. 101 070401
[25] Matsuoka T Shimizu K 2009 Nature 458 186
[26] Bazhirov T Noffsinger J Cohen M L 2010 Phys. Rev. 82 184509
[27] Marqués M McMahon M I Gregoryanz E Hanfland M Guillaume C L Pickard C J Ackland G J Nelmes R J 2011 Phys. Rev. Lett. 106 095502
[28] Matsuoka T Sakata M Nakamoto Y Takahama K Ichimaru K Mukai K Ohta K Hirao N Ohishi Y Shimizu K 2014 Phys. Rev. 89 144103
[29] Schaeffer A M Temple S R Bishop J K Deemyad S 2015 Proc. Nat. Acad. Sci. USA 112 60
[30] Saiz E G Gregori G Gericke D O Vorberger J Barbrel B Clarke R J Freeman R R Glenzer S H Khattak F Y Koenig M Landen O L Neely D Neumayer P Notley M M Pelka A Price D Roth M Schollmeier M Spindloe C Weber R L Woerkom L V Wünsch K Riley D 2008 Nat. Phys. 4 940
[31] Kugland N L Gregori G Bandyopadhyay S Brenner C M Brown C R D Constantin C Glenzer S H Khattak F Y Kritcher A L Niemann C Otten A Pasley J Pelka A Roth M Spindloe C Riley D 2009 Phys. Rev. 80 066406
[32] Bryk T Klevets I Ruocco G Scopigno T Seitsonen A 2014 Phys. Rev. 90 014202
[33] Chen M Vella J R Panagiotopoulos A Z Debenedetti P G Stillinger F H Carter E A 2015 AICHE J. 61 2841
[34] Vella J R Stillinger F H Panagiotopoulos A Z Debenedetti P G 2015 J. Phys. Chem. 119 8960
[35] Bakanova A A Dudoladov I P Trunin R F 1965 Sov. Phys.-Sol. State 7 1307
[36] Young D A Ross M 1984 Phys. Rev. B 29 682
[37] Su J T Goddard W A III 2007 Phys. Rev. Lett. 99 185003
[38] Kim H Su J T Goddard W A III 2011 Proc. Nat. Acad. Sci. USA 108 15101
[39] Rice M H 1965 J. Phys. Chem. Solids 26 483
[40] Thiel M V 1977 Compendium of shock wave data Lawrence Livermore Laboratory Report UCRL-50108 323
[41] Marsh S P 1980 LASL Shock Hugoniot Data Berkeley University of California Press
[42] Kresse G Hafner J 1993 Phys. Rev. 47 558
[43] Kresse G Furthmüller J 1996 Phys. Rev. 54 11169
[44] Nosé S 1984 J. Chem. Phys. 81 511
[45] Mermin N D 1965 Phys. Rev. 137 A1441
[46] Perdew J P Burke K Ernzerhof M 1996 Phys. Rev. Lett. 77 3865
[47] Blöchl P E 1990 Phys. Rev. 41 5414
[48] Anderson P W 1961 Phys. Rev. 124 41
[49] Gonze X Beuken J M Caracas R Detraux F Fuchs M Rignanese G M Sindic L Verstraete M Zérah G Jollet F et al. 2002 Comput. Mater. Sci. 25 478
[50] Mazevet S Zérah G 2008 Phys. Rev. Lett. 101 155001
[51] Recoules V Mazevet S 2009 Phys. Rev. 80 064110
[52] Holzwarth N A W Tackett A R Mattews G E 2001 Comput. Phys. Commun. 135 329
[53] Fabbris G Lim J Veiga L S I Haskel D Schilling J S 2015 Phys. Rev. 91 085111
[54] Kunz C Petersen H Lynch D W 1974 Phys. Rev. Lett. 33 1556